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Current  and  future  space-based  systems  have  rather  complex  mission  requirements  which  demand  the 
storage  of  large  amounts  of  liquid  propellants  on  board.  With  large  controller  bandwidths  and  rapid 
maneuvering  of  the  spacecraft  in  a  low  gravity  environment,  potential  coupling  between  the  sloshing  liquid, 
the  spacecraft  motion  and  structural  modes  need  to  be  carefully  evaluated  to  ensure  the  system  design 
adequacy.  For  achieving  the  mission  success,  the  first  important  step  is  to  understand  the  nonlinear  dynamics 
of  the  liquid  sloshing. 

The  repon  sununai  izes  a  two-year  study  on  the  development  and  application  of  the  numerical  method 
for  three-dimensional  liquid  sloshing  simulation.  Fluid  dynamics  and  fluid  loading,  including  total  force  and 
impact  for  the  ves.sel  undergoing  rapid  movement  were  simulated.  Effects  of  baffles  and  active  baffles  with 
or  without  feedback  mechanism  for  sloshing  control  were  compared.  It  was  found  that  moving  baffles  can 
be  very  effective  in  suppressing  large  amplitude  sloshing.  Complicated  swirling  intensification  by  drainage 
was  also  numerically  simulated. 
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NONLINEAR  SLOSHING  AND  THE  COUPLED  DYNAMICS 
OF  LIQUID  PROPELLANTS  AND  SPACECRAFT 


ABSTRACT 

Current  and  future  space-based  systems  have  rather  complex  mission 
requirements  which  demand  the  storage  of  large  amounts  of  liquid  propellants 
on  board.  With  large  controller  bandwidths  and  rapid  maneuvering  of  the 
spacecraft  in  a  low  gravity  environment,  potential  coupling  between  the 
sloshing  liquid,  the  spacecraft  motion  and  structural  modes  need  to  be 
carefully  evaluated  to  ensure  the  system  design  adequacy.  For  achieving  the 
mission  success,  the  first  important  step  is  to  understand  the  nonlinear 
dynamics  of  the  liquid  sloshing. 

The  report  summarizes  a  tV7o-year  study  on  the  development  and  application 
of  the  numerical  method  for  three-dimensional  liquid  sloshing  simulation. 
Fluid  dynamics  and  fluid  loading,  including  total  force  and  impact  for  the 
vessel  undergoing  rapid  movement  were  simulated.  Effects  of  baffles  and 
active  baffles  with  or  without  feedback  mechanism  for  sloshing  control  were 
compared.  It  was  found  that  moving  baffles  can  be  very  effective  in 
suppressing  large  amplitude  sloshing.  Complicated  swirling  intensification  by 
drainage  was  also  numerically  simulated. 
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CHAPTER  Is  INTRODUCTION 


Liquid  sloshing  in  a  moving  container  constitutes  a  broad  class  of 
problems  of  great  practical  importance  in  many  engineering  applications,  such 
as  tank  truck  accidents,  tank  car  derailing,  helicopter  flip-over,  liquid  fuel 
sloshing  in  space  rockets,  movement  of  liquid  cargo  in  ocean-going  vessels, 
and  on  a  much  larger  scale,  the  oscillations  of  water  in  lakes  and  harbors 
occurring  as  the  result  of  earthquakes.  The  potential  for  significant  loads 
resulting  from  sloshing  liquids  in  containers  has  been  realized  for  years. 
Such  loads  may  cause  structural  damage,  create  a  destabilization  effect  and 
produce  envirorimental  hazards  . 

Within  the  moving  container,  various  types  of  sloshing  waves  can  be 
created  depending  on  the  liquid  depth  and  motion  frequency.  In  an  unbaffled 
tank,  typical  forms  are  standing  waves,  travelling  waves,  solitary  waves, 
hydraulic  jumps  and  a  combination  of  these.  Associated  with  these  waves, 
impulsive  and  nori-impulsive  pressures  can  be  created  by  a  sloshing  liquid, 
Impulsive  pressures  which  are  usually  associated  with  hydraulic  jump  and 
travelling  waves  are  localized,  very  high  pressures,  with  a  duration  of 
1-10  msec,  caused  by  the  impact  of  the  liquid  on  the  solid  boundaries.  Non- 
impulsive  pressures  are  slowly  varying  dynamic  pressures  resulting  from 
standing  waves.  Usually  the  most  severe  iiispact  pressures  occur  on  the 
boundaries  near  the  stationary  liquid  level  or  at  the  corners  of  the 
container.  In  spite  of  harmonic  excitation,  these  pressures  are  neither 
harmonic  nor  periodic.  Various  baffles  produce  complicated  wave  forms,  most 
of  which  are  not  well  understood. 

The  problem  of  sloshing  has  been  studied  by  many  investigators  and  rather 
compreliensive  reference  lists  have  been  published  on  several  occasions^  ‘  ^ . 
Most  of  these  earlier  investigations  are  concerned  with  fuel  tanks  for  rockets 
for  the  space  program.  Much  of  the  sloshing  technology  developed  for  space 
applications  is  limited  to  spherical  and  cylindrical  containers,  v/hile  the 
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nature  of  slosh  loading  in  prismatic  tanks  is  probably  less  understood.  The 
emphasis  was  placed  on  frequencies  and  total  forces  as  they  relate  to  control 
system  requirements  and,  therefore,  the  effects  of  local  peak  impact  pressure 
on  structural  requirements  were  not  studied  to  any  extent.  Further,  the 
excitation  amplitudes  considered  in  space  applications  are  too  small  for  many 
other  simulations  such  as  ship  motion^.  Analytical  techniques  for  predicting 
large  amplitude  sloshing  are  still  not  fully  developed.  The  problem  is 
essentially  nonlinear  and,  therefore,  few  theories  are  available  for 
predicting  damping^ ' .  Our  principal  knowledge  of  damping  characteristic 
rema'j.ns  the  result  of  extensive  experimental  studies.  A  detailed  description 
of  sloshing  and  evaluation  of  slosh  loads  in  marine  applications  has  been 
given  in  a  report  by  Cox,  Bowles  and  Bass^'^.  Several  articles  of  fluid- 
svi  i'.iv.re  vibration  and  liquid  sloshing  in  reactor  technology  applications 
unciei  siav'thquake  excitation  can  be  found  in  Ma  and  Su^^ .  Sloshing  in  the  dam- 
rerervoi.r  sy.stem,  including  fluid- structure  interaction,  has  been 
Inves  tigated^^ ' . 

A  linear  mathematical  model  has  been  developed  for  rectangular  tanks 
executing  a  roll  oscillation  by  Graham^^  and  Chu,  et  al.^'.  In  their 
analytical  and  experimental  study  for  ship-roll  stabilization  tanks,  Chu,  et 
al.,  have  concluded  that  the  antirolling  tank  is  a  nonlinear  control  element 
throughout  its  practical  range  of  operation  and  that  a  nonlinear  mathematical 
model  must  be  developed  before  any  signific;ant  gain  over  present  design 
methods  can  be  foreseen. 

Nonlinear  theories  of  forced  oscillations  of  liquid  in  a  rectangular 
container  has  been  developed  by  Faltinsen^^  and  Verhagen  and  Wijiigaarden^^ . 
Both  studies  are  concerned  with  the  frequency  range  near  the  lowest  resonance 
frequency  and,  although  the  nonlinear  free  surface  condition  has  been  used, 
only  the  small  amplitude  roll  riir''ions  of  the  container  are  considered. 
Furthermore,  the  analysis  of  Verhagen  and  Wijngaarden  is  based  on  shallow- 
water  theory  and  its  application  is  limited  to  the  low  fill  depth  case,  while 
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the  analysis  of  Faltinsen  considers  the  depth  of  fluid  is  either  0(1)  or 
infinite.  The  comparison  between  theory  and  experiment  has  been  made  for  each 
study.  Both  authors  attributed  the  discrepancy  between  theory  and  experiment 
to  the  neglect  of  viscosity  iti  the.  theory.  Faltinsen  mentioned  that  "obvious 
noulinearit.ies  are  occurring  and  ic  Is  possible  chat  viscous  effects  play  a 
dominant  role.  .  .  further  study  must  investigate  Che  possibility  of 
incorporating  viscous  effects  in  an  approximate  way  in  our  potential  theory." 
In  a  later  paper^*^,  Faltinsen  developed  a  numerical  nonlinear  method  for  the 
study  of  sl.oshing  based  on  the  boundary  integral  technique.  For  small 
amplitude  liquid  sloshing  study,  the  finite  element  method  has  been  applied  by 
Ikekawa^-^  and  the  boundary  element  method  has  been  used  by  Nakayama  and 

Washizu22 . 

For  the  study  of  large  amplitude  sloshing  when  the  excitation  is  0(1)  , 
one  has  to  use  a  numerical  technique.  Numerical  techniques  have  been  used  to 
solve  time  -  dependent  incompressible  fluid  flow  problems  for  more  than  20 
years.  One  of  the  best  knowii  techniques,  the  Marker-and-Cell  (MAC)  method, 
uses  an  Eulerian  finite-di.fference  formulation  in  which  pressure  and  velocity 
are  the  primary  dependent  variables.  Hirt,  et  al.  of  Los  Alamos  Scientific 
Laboratory  developed  a  numerical  solution  algorithm  (SOLA) ,  using  a  finite 
difference  technique  based  on  the  MAC  method  for  solving  the  Navier-Stokes 
equations  for  an  incompressible  fluid'^^.  An  extension  of  the  SOLA  code,  SOLA- 
SURF,  permits  a  free  surface  to  be  located  across  the  top  of  the  fluid 
regions.  Navickas ,  et  al . ,  modified  the  SOLA-SURF  code  to  study  sloshing  of 
fluids  at  high- fill  levels  in  closed  tanks  to  predict  cargo  response 
characteristics  in  liquefied  natural  gas  tankers  at  high  loading  level.s  due  to 
both  periodic  and  raiidom  excitations^^ .  The  recent  sloshing  R  &  D  project  at 
Lloyd's  Register  of  Shipping  also  emphasizes  the  computer  analysis  of  sloshing 
I  effects  based  on  this  SOLA-SURF  code^^.  The  principal  restriction  of  the 

SOLA-SURF  code  is  that  the  free  surfaces  must  be  definable  by  s ingle -valued 
I  function.  Also,  the  slope  of  the  surface  must  not  exceed  Che  cell  aspect 

i 
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ratio  ay/ax.  In  marine  applications,  large  amplitude  excitations  are  often 
anticipated  with  violent  fluid  response  inevitable  and  these  limitations 
become  too  restrictive.  To  overcome  the  difficulty  of  tracking  the 
complicated  fluid  interface,  further  extension,  SOLA-VOF,  uses  the  volume  of 
fluid  technique26  to  track  a  free  fluid  surface. 

In  the  VOF  method,  the  volume  of  fluid  function  F  is  introduced,  defined 
as  unity  at  any  point  occupied  by  fluid,  and  zero  elsewhere.  The  discrete 
value  of  F  in  a  grid  cell  represents  the  fractional  volume  of  the  cell 

occupied  by  fluid.  The  distribution  of  F  contains  implicit  information  about 
the  location  of  the  iriterface  without  an  excessive  use  of  computer  time.  "Che 

multiple  interacting  free  surfaces  that  occur  in  liquid  sloshing  in  htfiled 

tank  are  simply  defined  by  the  VOF  method.  Su,  et  al.^^>^®,  using  the  moving 
coordinate  system  with  this  volume  of  fluid  technique,  studied  the  lU'i'.linear 
behavior  and  damping  characteristics  of  liquid  sloshing  in  partially  fil led 
enclosed  prismatic  tanks  subjected  to  a  large  amplitude  exc.ltation  lo''  mar  i  .le 
and  land  transportation  applications.  The  liquid  au.s  assumed  ':o  ne 
homogeneous  and  to  remain  laminar.  Sloshing  inside  [)art;.ally  f.i,]  '  d, 

enclosed,  baffled  tanks  and  unbaffled  tank.s  was  studied.  Se/eral  b..'.ff:e 
configurations  were  investigated  and  their  effects  on  sloshrng  have  keen 
identified.  A  few  comparisons  between  numerical  results  and  experimental  data 
have  also  been  made  with  good  agreement  observed  for  both  impact  and  nev 
impact  type  slosh  loads.  For  large  amplitude  excitations,  however,  the  liquid 
responded  violently  which,  after  an  initial  period,  caused  the  numerical 
solution  to  become  unstable.  In  an  ensuing  investigation,  an  improved  donor- 
acceptor  method  which  takes  into  account  surface  orientation  and  transports 
trapezoidal  shapes  from  cell  to  cell  is  used.  The  improved  algorithm  extends 
the  applicability  of  numerical  sloshing  simulation  to  permit  repeated  liquid 
impacts  on  the  surface  of  the  container ^ ^ .  Sufficient  number  of 
statistically  distributed  impact  pressure  can  then  be  generated.  Experimental 
investigation  was  carried  out  which  verified  the  accuracy  of  wave  height  and 
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Impact  ptessure  computation^^.  Very  recently,  a  three-dimensional  finite- 
difference  scheme  based  on  the  volume  of  fluid  technique  Vias  been 
developed^ ^ ^ ^ .  3-D  numerical  simulation  of  sloshing  in  arbitrary 
containers  becomes  feasible. 

Meanwhile,  current  and  future  spacecraft  mission  requirements  demand  the 
storage  of  large  amounts  of  liquid  propellants  on  board.  Space  applications 
again  attract  the  attention  of  sloshing  researchers ^ ^ .  Potential 
coupling  between  liquid  and  the  remaining  portion  of  the  spacecraft  can  be 
expected  due  to  large  liquid  mass  fr.acCions,  large  controller  bandwidths ,  low 
frequency  spacecraft  modes,  low-g  conditions  and  rapid  maneuverability 
requirements.  Understanding  the  sloshing  and  the  associated  phenomena  is 
essential  to  ensure  th.e  system  des'gu  adequacy  and  Che  mission  success.  The 
objective  of  the  research  was  to  extend  the  recent  progress  of  computational 
technique  to  understand  and  predict  the  dynamic  respon.se  of  liquid  in 
partially  filled  moving  containers  relevant  Co  spacccraf t/space  station 
applications  and  to  apply  this  technique  to  achieve  effective  fuel  tank  design 
optimization  and  baffle  selection. 
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CHAPTER  2.  VOLUME  OF  FLUID  TECHNIQUE  FOR  SLOSHING  SIMULATION 

The  case  to  be  considered  is  that  of  a  three-dimensional,  rigid,  enclosed 
tank  that  is  partially  filled  with  a  viscous  liquid.  The  liquid  is  assvmied  to 
be  homogeneous  and  to  remain  laminar.  A  general  motion  of  a  vehicle  -  fixed 
coordinate  system  (x,y,z)  with  respect  to  an  inertial  coordinate  system  is 
prescribed.  The  continuity  equation  for  an  incompressible  fluid  is  given  by 


—  (Bu)  +  ~(ev)  +  — (ew)  -  0 


The  equations  of  motion  are  the  Navier-stokes  equations  for  Newtoniati 
fluids 
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gy  —  body  acceleration  in  the  y  direction 
-  body  acceleration  in  the  z  direction 
The  gravitational  acceleration  is  included  as  part  of  body  accelerations, 
which  is  the  external  forcing  term  of  sloshing  motion. 

Finite  difference  representation  of  the  governing  partial  differential 
equations  is  obtained  as  described  in  Reference  [35],  The  finite  difference 
cells,  which  are  composed  of  cuboid  of  length  (5xj_,  width  <5yj  and  height  Sz^r, 
are  used  for  solving  the  governing  equations  numerically.  A  typical 
computational  grid  is  illustrated  in  Fig.  1.  The  tank  inside  is  divided  into 
the  IBAR  cells  in  the  x  direction,  JBAR  cells  in  the  y  direction  and  KBAR 
cells  in  the  z  direction.  A  single  layer  of  boundary  cells  surrounds  the  tank 
such  that  the  total  number  of  cells  is  IIIAX  (“IBAR+2)  times  JMAX  (“JBAR+.7) 
times  KMAX  (“■KBAR+2).  These  fictitious  cells  are  used  to  set  regular  boundary 
conditions  so  that  the  same  difference  equations  used  in  interior  cell  can 
also  be  used  at  the  boundaries.  A  staggered  grid  is  used  as  shov;n  iti  Fig.  2 
with  u- -velocity  at  the  right  face  of  the  cell,  w- -velocity  at  the  front  face 
of  the  cell,  v- -velocity  at  the  cop  face  of  the  cell,  and  pressure  and  F  '  he 
vol>ime  of  fluid  function)  at  the  cell  center.  Fig.  2  also  shows  the  locations 
of  the  geometrical  quantities  AR,  AT,  AF.  and  VC,  arising  from  the  partial 
cell  treatment,  which  are  the  fractions  of  the  right  cell  face,  top  cell  face, 
front  cell  face  and  cell  volume  open  to  flow.  Geometrical  arrays  AR ,  AT,  AF, 
and  VC  must  be  defined  if  an  obstacle  or  a  curved  boundary'  is  included  in  the 
me.sh . 

The  finite  difference  approximation  representing  the  continuity  equation 
for  a  typical  cell  (i,j,k)  is: 
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The  finite  difference  approximation  of  momentum  equations  are; 
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Where  the  convective  and  viscuous  fluxes  in  x  direction  are  defined  as; 
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^i  ^i+1 


(5x.  ,  (5x. 

,  n  ,<■  s  ,  /  n  ,  1+1  ,  n  n  .  i 


6x . 

1 


6x.  . 
1+1 


-  6X^  f  k  <.*xgn(u^  J  J.)  .  «x.! 


/  *,r  s  i+1/2  ,  n  n  .  *^^j-l/2  ,  n  n  ^ 

'  ^u"  ,  i.,J,k  i,j-l,k  i,j+l,k  i,j,k 

.^j-1/2  ^j+1/2  J 


/  *,c  N  .  /  *^^j+l/2  /  n  n  *^^j-l/2  n  n  ^  /inx 

(V  /Jyjax.s„(v  )  - (u.  .  i^-u.  ,.J  1^)  -  - -  (“i ,  j,l ,  k"*!,  J  ,  k'  O") 


^U/2  - 


-  ‘yj+1/2  ^  '’h-1/2  *  -  ■Syj.i/j) 


1  r  1  X  '1  ri  .  ^  ,  n  ,  n  . 

- dx.(v.  ,  .  ,+v.  ,  .  ,  ,  )  +  i5x.,t(v.  .  +v.  .  ) 

1  l-xl.J.V.  1+1  l.J.k  l.J-l.k 

1  1  +  1 


2  ((5x . 


r-,.-x  .  S  ‘^^k+1/2  ,  n  n  ,  ^^k-1/2  ,  n  n  ,  ^ 

FIjZ  “  (w  /6z  ) - - —  (u.  .  -u.  1  )  +  -  (u.  .  ,  . ,  -u.  )  + 

^  u/  i,j,k  i,j.k-l'  i,J,k+l  i,j,k 

.  ^kl/2  ^^k+l/2  J 
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(w  /(5z^)a*sgn(w  ) 


(5z 


6z 


k+1/2  ,  n  n  ,  '^^k-1/2  ,  n  n 

- - —  (u.  .  ,  -u.  - -  (u.  .  ,  n  -U,  .  , 

l.J.k  l.J.K-1  l,J.k+l  l.J.k 


k-1/2 


k+l/2 


‘"k-1/2  •  ’  ^"k-l>''^ 


‘"u  -  *\n/2  ■"  '*'k-l/2  +  »*=8n<"  >  <‘^k+l/2  '  '*"k-l/2> 


w 


2(6x.+5x.^J 


_.ri  i*i  ve/^i  T*i  V 

OX.(w.  ^  .  ,+W.  -  .  ,  ,  )  +  ox.  t(w.  .  ,+W.  .  ,  t) 

1  1+1, j.k  i+l,j ,k-l  1+1  i.J.k  i,j,k-l 


VISX  -  V  (  +  —  +  -—  ) 

ay'-  az^ 


;)2 

a  u 


.2  ax.  ,  ,oax. , 1  ax. 

ax  1+1/2  1+1  1 


.  ,  n  n  ,  .  n  n  , 

ax. (u.  ,  .  ,  -u.  .  , )  -  ax.  , (u.  .  ,  -u.  ,  .  ,  ) 

1^  1+1, j,k  i.j.k  1+1  1,1, X  i-l,j,k 


a  u 


^yj+l/2^'^^j+l/2  ^^j-1/2 


,  n  n  .  ,  ,  n  n  . 

6y .  .  ,„(u.  .  ,  -u.  .  ,  )  -  ay.  T  .  ,  -u.  .  ,  ,  ) 

-’^j-i/2^  i,j+l,k  i,j,k  ■^J+1/2  i,J,k  i,j-l,k 


a^u 


az^  '^^k-l/2^  ‘’^k+l/2^'^^k4l/2 


c  /II  n  \  j;  /  n  I* 

oz,  T  ,,, (u.  .  ,  .  -  VI.  .  ,  ;  -  az,  .  ._(u.  .  -u.  .  ,  T 

k-1/2^  i,j,k(l  i,j,k  kH 1/2  i.j.k  i,j,k-l 


In  the  equal  mesh  condition, 


ax.  T  -  ax.  a.x.  ,  -  ax 

1+1  1  1-1 


«yj.,i  -  Oyj  -  ■Syj.i  -  «y 


az,  .  -  az,  "  az.  ,  -  az 

k+ 1  k  k - 1 


)  (11) 


(12) 


(13) 


(14) 


(15) 
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The  expressions  for  FXJX,  FXJY,  FUZ,  VISX  can  be  simplified  as  follows: 


FUX 


n 

_  j 

26k 


n 


A  -  a*sgn(u  ) 

» J  I  ^ 


FlTf  - 


V 


2dx 


B  “  a*sgn(v  ) 


1  ,  n  n  n  n  , 

- (v,  ,  +  V.  ,  .  T  ,  +  V.  .  ,  +  V.  .  ,  ,  ) 

1+1, j,k  i+l,j-l,k  i,J,K  i,j-l,k 


FUZ  - 


26x 


^'■^>"i.J,k+l  "  K.j.k  -  (l"">^.j.k-ll 


C  -  a*sgn(w  ) 
+ 

w  -  - 

4 


1  ,  n  n  n  n  . 

—  (w. , ,  .  ,  +  w.  T  ,  ,  +  w.  .  ,  +  w.  .  ,  T ) 
1+1, j,k  i+l,j,k-l  i,j,k  i,j,k-l 


(16) 


(17) 


(18) 


VISX  -  -  (u.  T  -  2u.  .  ,  +  u.  T  ) 

„  1+1 ,J ,k  i,J,k  1-1, j,k 

(5x 


u  ,  u  „  n  n  . 

+  -  (u.  .  ,  ,  -  2u.  .  ,  +  u.  .  .,  ,  ) 

„  i,J+l,k  i,j,k  i,j-i,k 


6y 


u  ,  XI  „  n  n  . 

+  -  (u,  .  ,  ^  -  2u.  ,  ,  +  u.  .  ,  ,/ 

1  ,  j ,k+l  1 , J ,k  1 , j ,k-l 

6z 


(19) 


As  the  same  way,  we  can  define  the  expressions  for  convective  and  viscous 
fluxes  in  y  and  z  direction,  FVX,  FVY,  FVZ,  VISY,  FWX,  FWY,  FWZ  and  VISZ. 

The  coefficient  a  in  these  expressions  gives  the  desired  amount  of 
upstream  (donor  cell)  differencing;  that  is,  when  a  is  zero,  these  difference 
equations  are  second  order  accurate,  centered  difference  approximation.  Wlien 
Q  is  equal  to  unity,  the  equations  redvice  t.o  the  full  upstream  or  donor  cell 
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form.  The  proper  choice  for  a  is  needed  to  in.sure  numerical  stability  and  is 
discussed  later. 

In  order  to  satisfy  these  equations,  the  well-known  SOIA-VOF  algorithm  is 
employed  to  satisfy  the  three  difference  equations  in  each  control  volume  and 
to  correct  the  results.  At  each  time  step,  we  must  solve  for  the  u- -velocity. 
V- -velocity,  w- -velocity  and  pressure  in  each  control  volume.  Explicit 
approximations  of  the  velocity  field  are  obtained  from  the  momentum  equations, 
Eq.  (6), (7), (8),  using  old  time  level  values  for  the  advective,  pressure, 
viscous,  and  body  force  terms.  In  order  to  satisfy  the  continuity  equations, 
the  pressure  (and  velocities)  must  be  adjusted  in  each  computational  cell 
occupied  by  fluid. 

The  finite-difference  form  of  continuity  equation  (5)  can  be  rewritten  as 


S  - 


VC.  .  . 
i.J  .k 


1  ,AT.  .  ,  i,)  +  -^  V  lAF.  .  , ) 

i,j-l,k  r,j-l,k  .  1.1, k  i.j.k  i,j,k-l  i,j,k-l 


(20) 


This  equation  is  an  implicit  relation  for  the  new  pressures.  A  solution  may 
be  obtained  by  the  iterative  process.  In  each  cell  containing  fluid,  but  not 
a  free  surface,  the  pressure  change  needed  to  drive  the  S  of  Eq .  (20)  toward 
zero  is 


(5p  =  _1?_  (21) 

as/dp 

where  S  is  evaluated  v/ith  the  most  updated  values  of  p,u,v,  and  w  available, 
and  the  derivative  is  with  resjiect  to  Pi^j^l^-  The  equation  (21)  is  derived  by 
substituting  the  Eqs .  (25)  through  (30)  below  into  equation  (20)  and  solving 

for  i5p.  The  quantity  /3  =  l/(as/ap)  is 


VC 


i ,  j  ,  k 


2it(A.+A.  _^+^.+f 


(22) 
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whex’e 


X. 

1 


AR.  .  , 

i.J  ,k 


<5x.  (6x.+5x.  , ) 
1  1  1+1 


AR. 


i-1.  j  ,k 


<5x.  ((^x.+5x.  , ) 
1  1  1-1 


AT.  .  , 

i.J 


^  <57^  ('5yj+'5yj^.i^) 


^j-1 


AT.  .  ,  , 

1  ■  J  - 1  ■  V<'- 

(5y.  ((5y  +(5y  ) 

J  J  J 


AF.  .  , 
_ i. J ,k 


and 


^k-1 


AF.  ... 

i.J  .tc-1 


In  equal  mesh  condition 


VC.  .  , 
i,J  .k 


6t 


^^i,  j  ,k‘^^i-l,  j  ^^^i,  j  .k"^^^!,  j J'^^i  ,  j  .k'^^^i  ,k.k-l 

^x  6y  6z 


)  (23) 


The  new  estimate  for  the  cell  pressure  is  then 


n.j.k  ^ 


(24) 


and  new  estimates  for  the  velocities  located  on  the  sides  of  the  cell  are 


u .  .  ,  + 

i.J  ,k 


<5t(5p 
(5x . 


(25) 


u. 


i-1 .  j  ,k 


(5t6p 

(5x . 

1 


(26) 
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V.  + 

I.J  .k 


<5t(5p 

^yj 


j -l.k 


(5t6p 

<5yj 


w.  .  ^ 

<5z^ 

k 


w.  + 

i.J  ,k-L 


6t(!)p 

6z, 


(27) 


(28) 


(29) 


(30) 


This  procedure  is  modified  for  cells  containing  a  free  surface.  For 
these  cells,  the  free  surface  boundary  condition  is  satisfied  by  setting  the 
surface  cell  pressure  Pi  j  ^  equal  to  the  value  obtained  by  a  linear  inter¬ 
polation  between  the  pressure  wanted  at  the  surface  p^  and  a  pressure  of  its 
interpolation  cell  inside  the  fluid  so  the  S  function  is: 

S  -  (l-f?)Pn  +  »7Ps  ■  Pi,  j  ,k 

where  r]  -  d^/dj,  is  the  ratio  of  the  distance  between  the  cell  centers  and  the 
distance  between  the  free  surface  and  the  center  of  the  interpolation  cell  as 
shown  in  Fig.  3.  A  (5p  is  obtained  from  the  new  pj^  and  the  old  Pi,j_i^  and  the 
new  j  is  obtained  interactively  by  under- relaxation. 

Before  using  Eq .  (31)  for  surface  cell  interpolation,  we  must  first 
determine  the  distance  dg ,  the  distance  between  the  free  surface  and  the 
center  of  the  interpolation  cell.  Since  the  free  surface  is  assumed  either 
nearly  horizontal  or  vertical  in  original  VOF  method,  the  real  surface  in 
Fig.  4(a)  is  treated  as  the  surface  in  Fig.  4(b).  Therefore,  the  free  surface 
condition  does  not  impose  the  right  position  and  large  error  can  occur  due  to 
this  approximation.  TViis  situation  may  become  even  worse  in  three-dimensional 
problem,  the  free  surface  of  which  is  far  more  complex.  The  algorithm  of 
determining  distance  dg  is  improved  in  Reference  [35].  The  free  surface  is 
treated  as  a  plane.  The  slopes  of  plane  AN  and  BN  are  computed  by  the  F  value 
of  its  adjacent  cells  and  discussed  in  the  next  section.  Tlie  distance  dg  is 
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determined  by  letting  the  volume  between  the  plane  and  bottom  face  be  equal  to 
the  F  value  of  surface  cell.  Fig.  5  illustrates  the  position.^  of  two  kinds  of 
approximate  plane  in  surface  cell.  In  the  case  of  Fig.  5(a),  the  plane  is 
entirely  in  the  surface  cell  and  F  value  in  surface  cell  satisfies 


F  > 


(32) 


where,  6x  -  (5x/(5y  and  Sz  -  6z/8y .  The  distance  d^  is 


d  -  (0.5  +  F  ,  )6y  (33) 

^  ^  I  J  1 

In  the  case  of  Fig.  5(b),  the  approximate  plane  intersects  the  bottom  face  of 
surface  cell  and  F  value  in  surface  cell  satisfies 


^  BN—  AN-r- 
F  <  — (5x  +  — 6z 


(34) 


Through  geometrical  analysis  in  this  case,  we  can  obtain  the  distance 
follows ; 


d  - 

s 


^  BN^—  ,  ANv-,  7 

(0.5  +  — (5x  +  — (5z)  -  <5 


(35) 


<5  in  Eq.  (35)  is  one  of  the  roots  of  the  following  cubic  algebraic  equation 


-3  - 

6  +  pd  +  q  -  0 


(36) 


where 


p  -  -6(AJ^)  (BN)  «5x  6z 


_  _  DM _  AM _ 

q  -  6(AN)  (BN)  (5x  6z  (— (5x  +  —dz  -  F.  .  ) 

O  O  t  .  J  > 


A  complete  iteration,  therefore,  consists  of  adjusting  pressures  and 
velocities  in  all  cells  occupied  by  fluid  according  to  Eq .  (24)  through 
Eq.  (30).  Convergence  of  the  iteration  is  achieved  when  all  cells  have  S 
values  whose  magnitudes  are  below  some  small  number,  c.  Typically,  e  is  of 
order  10’^,  although  it  can  vary  witli  the  problem  being  solved  and  the  units 
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chosen  for  the  problem.  The  last  iterated  quantities  of  velocity  and  pressure 
are  taken  as  the  advanced  time  values . 

In  some  cases,  convergence  of  iteration  can  be  accelerated  by  multiplying 
<5p  from  Eq,  (21)  by  an  over- relaxation  factor  w.  A  value  of  w  -  1.7  has 
provided  relatively  efficient  results,  ’'Ut  the  optimal  value  is  flow- 
dependent.  A  value  of  two  or  greater  gives  an  unstable  iteration.  In 
practice,  the  free  surface  condition,  Eq.  (31)  leads  to  an  over-relaxation 
type  of  instability  when  the  interpolation  factor  w  is  greater  than  one. 
Stability  can  be  insured  by  under  -  relaxing  the  pressure  variations  in  cells 
used  as  interpolation  neighbors  for  surface  cells. 

Tracking  the  free  surface  is  essential  in  a  free  surface  flow 
computation.  The  free  surface  is  treated  by  introducing  a  function  F(x,y,z,t) 
that  is  defined  to  be  unity  at  any  point  occupied  by  the  fluid  and  zero 
elsewhere.  Thus,  F~1  implies  a  cell  full  of  fluid,  while  F-0  denotes  an  empty 
cell.  A  cell  with  F  values  between  zero  and  one  are  partially  filled  with 
fluid;  they  are  either  intersected  by  a  free  surface  or  contain  voids 
(bubbles)  smaller  than  cell  mesh  dimensions.  A  free  surface  cell  (i,j,k)  is 
defined  as  a  cell  containing  a  non-zero  value  of  F  and  having  at  least  one 
neighboring  cell  (i±l,j,k),  (i,j±l,k),  (i,j,k±l),  that  contains  a  zero  value 
of  F.  The  F  function  is  utilized  to  determine  which  cells  contain  a  boundary 
and  where  the  fluid  is  located  in  those  cells.  Additionally,  the  derivatives 
of  F  can  be  used  to  determine  the  mean  local  surface  normal,  and  using  also 
the  cell  F  value,  to  construct  a  plane  cutting  the  cell  that  will  approximate 
the  interface . 

The  time  dependence  of  F  is  governed  by 


SF  dF  dF  3F 

— -  +  U -  +  V —  +  w - 

3t  dx  dy  dz 


0 


(37) 


We  combine  Eq .  (37)  and  continuity  equation  (1)  to  obtain 

i-(eF)  +  — (euF)  +  — (evF)  +  --(ewF)  =  o 

5t  dx  dy  dz 


(38) 
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where  6  is  the  partial  cell  parameter  defined  as  before.  The  difference 
approximation  of  Eq .  (38)  is 


-  F 

i,j,k  i.j.k 


(5t 


VC.  .  , 
i.J  .k 


(5x. 

1 


(AR.  .  ,F.  .  -  -AR.  ,  .  ,  u"'*'?'  .  ,  F. 

i,j ,k  1.3 ,k  L.j ,k  i-l.j ,k  1-1,3 ,k  i-l.j ,k" 


(AT.  .  ,  F.  .  ,  -AT.  .  ,  ,  ,  F.  .  ,  ,  ) 

i,3,k  i,3,k  i.j.k  i.j-l.k  1,3-1, k  1,3-1, k 


(aF  F  -AF  F  ) 

^  i.3.1^  i.j.k  i,j,k  i,j,k-l  i.j.k-l  i.j.k-l'' 


(39) 


which  serves  as  the  basis  for  the  convection  of  F. 

The  convection  algorithm  must  (1)  preserve  the  sharp  definition  of  free 
boundaries;  (2)  avoid  negative  diffusion  truncation  errors;  and  (3)  not  flux 
more  fluid,  or  more  void,  across  a  computing  cell  interface  than  the  cell 
losing  the  flux  contains.  To  accomplish  this,  a  type  of  donor-acceptor  flux 
approximation  from  original  SOL.A-VOF  is  improved  and  employed  in 
Reference  (35).  In  this  FAU  modified  VOF  technique,  it  is  assumed  that  the 
fluid/void  interface  is  a  plane  segment  cutting  through  the  surface  cell.  So 
all  the  surface  physics  algorithms  must  be  improved. 

The  basic  orientation  of  a  surface  cell  is  described  by  the  parameter  NF 
through  its  numerical  values  1,  2,  3,  4,  5,  or  6.  The  definition  of  numerical 
values  for  the  NF  is  specified  in  Table  1.  The  differential  geometry 
describing  the  svirface  curvatures  depends  strongly  on  the  numerical  value 
selected  for  NF;  it  is  important  to  make  an  appropriate  choice  for  NF. 

Nevertheless  it  has  not  been  possible  to  devise  a  con.pletely 
satisfactory  algorithm  of  general  applicability  of  the  choice  of  NF  values. 
The  algorithms  that  are  currently  implemented  in  the  FAUMVOF  code  suffice  to 
provide  reasonable  simulations  for  the  problems  to  which  the  code  has  been 
applied  so  far. 
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TABLE  1 


DEFINITION  OF  VALUES  FOR  NF 


NF-0  fluid  cell,  contains  fluid  and  has  no  void  adjacent  to  any  of  its  faces. 
NF-1  surface  cell,  fluid  most  nearly  on  the  left  of  the  surface  cell. 

NF-2  surface  cell,  fluid  most  nearly  on  the  right  of  the  surface  cell. 

NF'-S  surface  cell,  fluid  most  nearly  under  the  bottom  of  the  surface  cell. 

NF=4  surface  cell,  fluid  most  nearly  above  Lhe  top  of  the  surface,  cell. 

NF“5  surface  cell,  fluid  most  nearly  in  the  back  of  the  surface  cell. 

NF“6  surface  cell,  fluid  most  nearly  in  the  front  of  the  surface  cell. 

NF--7  isolated  cell,  contains  fluid  but  all  cells  adjacent  to  one  of  its  faces 
are  void. 

NF=8  void  cell,  contains  no  fluid. 


We  calculate  the  value  of  NF  by  deciding  where  the  fluid  is  mostly 
located.  If  surface  cell  has  only  one  empty  adjacent  cell,  the  value  of  NF  is 
made  by  choosing  the  adjacent  cell  opposite  to  empty  cell,  that  cell  is  also 
the  interpolation  neighbor  fluid  cell.  For  example,  the  value  of  NF  is  3  and 
cell  (i,j-l,k)  is  the  interpolation  neighbor  fluid  cell,  if  cell  (i,j+l,k)  is 
only  empty  cell  of  surface  cell  (i,j,k).  If  surface  cell  has  two  empty 
adjacent  cells,  the  value  of  NF  is  made  by  choosing  from  the  two  adjacent 
cells  opposite  the  empty  cells,  that  cell  with  the  largest  F  value,  i.e., 
containing  the  most  fluid,  to  be  the  interpolation  neighbor  fluid  cell.  For 
example,  if  two  adjacent  cells,  cell  (i+l,j,k)  and  cell  (i,j+l,k)  are  the 
empty  cells  and  F  value  of  cell  (i-l,j,k)  is  greater  than  that  of  cell 
{i,j-l,k),  the  value  of  NF  is  1  and  cell  (i-l,j,k)  is  chosen  as  interpolation 
neighbor  fluid  cell. 
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The  two  slopes  of  the  interface  AN  and  BN  can  be  defined  according  to  the 
value  of  NF  and  are  discussed  at  the  end  of  this  section.  The  interface  is 
approximately  represented  by  a  plane  equation  at  its  local  coordinates  ^x^y^'z  ■ 
For  example,  if  value  of  NF  is  3,  the  local  coordinate  is  such  that  its  origin 
is  at  the  center  of  bottom  face  of  the  surface  cell  and  is  upward.  The 
plane  equation  is  given  by 


^y  '0  X  z 


two  restrictions  for  (  are 

y 


0  <  c  <  6y. 

y  J 


^  “  (5y . 

y  J 


^  -  0 
y 


^  >  <5y . 
y  J 


C  <  0 

y 


where  70  distance  between  the  free  surface  and  bottom  face  of  the 

surface  cell.  From  the  definition  of  r]  —  d^^/d^  in  Eq.  (31),  we  can  obtain  7q 


70  -  dc(l/r?  -  0.5)  •  (43) 

The  improved  donor-acceptor  method  is  employed  to  calculate  the  advection 
of  F  through  the  cells.  First,  at  each  face  of  each  computing  cell,  the  two 
cells  immediately  adjacent  to  the  interface  are  distinguished,  one  becoming  a 
donor  cell  and  the  other  an  acceptor  cell,  and  cell  quantities  are  given  the 
subscripts  D  and  A,  respectively,  e.g.,  Fj),  F^.  The  labeling  is  accomplished 
by  means  of  the  algebraic  sign  of  the  fluid  velocity  normal  to  the  face;  the 
donor  cell  is  always  upstream,  and  the  acceptor  cell  downstream,  of  the  face. 
We  emphasize  that  the  D  and  A  labels  are  assigned  separately  for  each  cell 
facz.  Thus,  each  computational  cell  will  have  six  separate  assignments  of  D 
or  A  corresponding  to  each  of  its  cell  face. 

The  amount  of  F  fluxed  across  the  donor  cell  face  in  one  time  step  can  be 
calculated  for  different  value  of  NF.  Fig.  6  gives  an  illustration  of 
F-advection  in  the  typical  case  when  NF=3.  In  this  case,  the  amount  of  F 
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across  the  right  face  of  the  cell  at  one  time  step  (5t  is  AFj^  shown  in 

Fig.  6(a).  By  the  partitioning  of  the  into  M  equal  subdivisions,  each  of 

length  5^2'  with  midpoints  ^2! . €zm . ^zM'  obtain  the  as 


AF 


U(5t 

6v .  6z, 

J  k 


Z  ,  ^  <5x.  --  u5t.^  )(5f 

ra-1  y  X  ^  ^2111 


(4A) 


The  amount  of  F  fluxed  across  the  top  face  of  the  cell  at  one  time  step  is  AFy 
shown  in  Fig.  6(b).  Again  by  the  partitioning  of  the  5xj  into  L  equal 

subdivisions,  each  of  length  and  the  midpoints  ^  . ^  xl  . ^xL' 

can  obtain  the  AFy  as 

- ^ - h  (45) 

y  j.  t  T-'m-l  y  ^x  z  '  ^ 

■'  6x.(5z,  •' 

1  k 


where 


C  -  MAX  k.(f  ,^_,_)-(6y,-V(St)  .  0 

>  1  >  Xi  ,1,11.  j 


(46) 


The  amount  of  F  across  the  front  face  of  the  cell  at  one  time  step  is  AF 
shown  in  Fig.  6(c).  We  can  obtain  the  AF^  as 


z 


AF  - 
z 


w6t 

<5x.  <5y . 
1  J 


(47) 


To  prevent  the  fluxing  of  more  fluid  from  the  donor  cell  than  it  contains, 
the  following  restrictions  are  applied  to  4F^,  AF  and  AF^ 


AF  -  MIN  1 

X  1 

,  j  ,  k] 

(48) 

AF  -  MIN  1 

y  1 

[“y 

^D'VS/^"i,j,k] 

(49) 

AF  -  MIN  1 
z  { 

V^D"V^"i,j.k] 

(50) 

The  complete  flu.xing  algorithm  is  applied 

computing  cell  faces.  When  the  necessary  fluxes 

independently  at  the  six 

have  been  computed,  F  is 

advanced  through  one  time  step  using  Eq .  (39). 
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Sometimes,  ;5bove  F-advection  algorithm  will  generate  spurious  small  wisps 
of  fluid  in  the  void  cells  of  the  computing  mesh.  An  algorithm  used  to 
suppress  these  spvirious  wisps  of  fluid  is  to  set  a  lower  bound  for  F, 
i.g.  0.05.  If  a  surface  cell  is  trying  to  flux  material  into  an  empty  cell, 
the  flux  is  sec  to  zero  until  it  is  greater  chan  0,05.  The  limiting  value  of 
0.05  may  not  be  optimal  for  all  problems. 

Truncation  errors  and  rounding  errors  can  cause  F-values  determined  by 
the  above  procedure  to  occasionally  have  values  slightly  less  than  zero  or 
slightly  greater  than  unity.  Therefore,  after  the  advection  calculation  has 
been  completed,  a  pass  is  made  through  the  mesh  to  reset  values  of  F  less  than 
zero  back  to  zero  and  values  of  F  greater  than  one  back  to  one.  Accvimulated 
error  in  fluid  volume  introduced  by  these  adjustments  are  recorded  and  may  be 
printed  at  any  time. 

There  is  a  final  adjustment  needed  in  F  so  chat  it  may  be  used  as  a 
boundary  cell  flag.  Boundary  cells  have  values  of  F  lying  between  zero  and 
one.  However,  in  a  numerical  solution,  F-values  cannot  be  tested  against  an 
exact  number  like  zero  atid  one  because  rounding  rounding  errors  would  cause 
spurious  results.  Instead,  a  cell  is  defined  to  be  empty  of  F  when  F  is  less 
than  EMF  and  full  when  F  is  greater  than  1-EMF,  where  EMF  is  typically  10"^. 
If,  after  advection,  a  cell  has  an  F  value  less  than  EMF,  this  F  is  set  to 
zero  and  all  neighboring  full  cells  become  surface  cells  by  having  their 
F-values  reduced  from  unity  by  an  amount  1.1*EMF.  These  changes  in  F  are  also 
included  in  the  accumulated  volume  error.  Volume  errors,  after  hundreds  of 
cycles,  are  typically  observed  to  be  a  fraction  of  1%  of  the  total  F  volume. 

Following  the  calculation  of  the  new  F-values  for  all  cells,  the  new  cell 
types  are  redefined  and  appropriate  flags  NF  are  assigned.  At  the  same  time, 
the  approximate  orientation  of  the  fluid  in  each  surface  cell  is  determined 
and  a  pressure  interpolation  neighbor  cell  is  assigned.  The  slope  of  inter¬ 
face  is  estimated  by  introducing  six  surface -he igb.t  functions  dY/SZ,  dY/dX, 
dZ/9X,  dZ/dY,  dX/dY,  dX/dZ  based  on  the  values  of  F  in  the  surface  cell  and 


its  neighbors.  The  good  approximations  to  these  functions  are: 


1 - wm  I  i.i  II  I  II  I 


\,k+l''^i.  j  -l.k+l'^-^j  -l'*'^i,  j  ,k+l^^j'^^i,  j-^l.k+l^^j+1 
’^i.k-l’^i,  j  -l,k-l‘*^j-l'^^i,  j  ,k-l^^j‘*'^i,  j+l.k-l'^^j  -1 


^  .k'^^i-l  .k^ 

ax  ax.  ,+2ax.+6x.  , 

1+1  1  1-1 


^  i+l ,  k  ^i+1 .  j  - 1 ,  k'^^j  -  l^^i+1 ,  j  ,  ^i+1 ,  j+1 .  k*^^j  +  l 

^i-l,k”^i-l.  j  -  l.k'^^j  -  l'^^i-1  .  j  ,k‘^-’'j'^^i-l,  j+l,k^^j+l 


ax  ^x.^,+2ax^+ax..^ 


2.  .  -“F..,  .  ,  .az,  ,+F.^,  .  ,az, +f.^t  .  az 

1+1, J  i+l,j,k-l  k-1  1+1, j,k  k  i+l,j,k+l  k+1 

z.  .-F.  ,  .  ,  .az,  ,+F.  T  .  ,a  ,+F.  1  .  ,  ,iaz,  . 

1-1, J  1-1, J, k-1  k-1  1-1, j,k  ic  i-l,j,k+l  k-J-l 


2CZ  -Z  ) 

^  >  •  i, j+1  i. j 

^i,  j+l’^i,j+l,k-l‘^^k-l^^i,j+l,k'^^k'‘'^i,j+l,k+l‘^^k+l 

^i,j  -1  ^i,j  -  l,k-l'^^k-l^^  i,j  -  l,k^^k^^i,j  -  l,k+l*^‘'k+l 


ax  ^^^j+l,k"^j-l,k^ 


ay  ay.  ,+2ay.+ay.  , 

■^J+1  J  ^J-1 


X.  T  ,-F.  ,  .  ,  ,  ax.  ,+F.  ,  ax.+F,  -  . ax.  , 

j+l,k  i-l,j+l,k  1-1  i,j+l,k  1  i+i,j+i,k  1+1 

X.  ,  ,“F.  ,  .  .  ,ax.  ,+F.  ,  ,  ,ax.+F.^,  .  ,  ,ax.^, 
j-l,k  1-I,ji.,k  1-1  i,j-l,k  1  1+I,j-I,k  1+1 


E  .  ^-"j,^+l'"j.k-l^ 

az  az,  ,  +2z,  +z,  , 
k+1  k  k-1 


(51) 


(52) 


(53) 


(54) 


[55) 


(56) 
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X.  ,  1-F,  ,  .  ,  ,i5x.  ,+F.  .  ,  ,5x.+F.  ,  .  ,  ,6x. 
j,k+l.  i-l,j  ,k+l  1-].  i,j,k+l  1  i+l,j.k+l  1+1 

X.  ,  t“F.  t  .  ,  ^<5x.  1+F.  .  (5x.+F.  T  .  ,  t6x.  _ 

j,k-l  i-l,j,k-l  1-1.  i,,i,k-l  1  i+l,j.k-l  1-1 

For  each  surface  having  certain  NF  value,  the  slope  of  the  interface  AN. and  BN 
can  be  assigned  according  to  the  definition  in  Table  2. 


TABLE  2 

THE  DEFINITION  OF  VALUES  FOR  AN  AND  BN 


NF  »  1,2 

AN  =  ^ 

BN  = 

ax 

oY 

az 

NF  »  3,4 

AN-^ 

BN  - 

av 

az 

ax 

NF  -  5,6 

AM 

AN  -  — 

BN  - 

az 

ax 

av 

In  addition  to  the  free  surface  boundary  condition,  it  is  necessary  to 
set  conditions  at  all  mesh  boundaries  and  at  surface  of  all  internal 
obstacles.  At  the  mesh  boundar  ie..s ,  a  variety  of  conditions  may  be  set  using 
the  layer  of  fictitious  cells  surrounding  the  mesh.  Consider,  for  example, 
the  left  boundary;  if  this  is  a  rigid- free  slope  wall,  the  normal  velocity 
there  must  be  zero  and  the  tangential  velocity  should  have  no  normal  gradient, 
i  .  e . 

-  0 

VI, j ,k  “  V2,j ,k 

wi,j ,k  ”  W2, j  ,k  (57) 

Pl,j,k  "  P2,j  ,k 
Fl,j ,k  “  t2,j  .k 
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If  the  left  boundary  is  a  no-slip  rigid  wall,  then  the  tangential 


velocity  component  at  the  wall  should  also  be  zero, 
i .  e . 

n,j ,k  “  -^2,j ,k 

n,j  ,k  -  -W2,j  ,k  (58) 

Pl,j,k“  P2.j,k 
Fl,j,k“  f2,j,k 

Boundary  conditions  similar  to  those  for  the  left  wall  are  used  at  the 
right,  front,  back,  top  and  bottom  boundaries  of  the  mesh.  In  the  case  of  a 
tank  with  baffles,  velocities  are  sec  to  zero  all  the  time  for  the  interior 
obstacles  and  the  other  variables  are  maintained  at  their  initial  values. 

In  practice,  non- rectangular  containers  such  as  cylindrical  and  spherical 
containers  are  coiraiionly  used.  The  boundary  treatment  for  rectangular  tank 
wall  will  encounter  severe  difficulties  when  curved  boundary  is  encountered. 
For  a  free -slip  condition  on  tank  wall,  the  following  must  occur:  (a)  the 
velocity  normal  to  a  boundary  surface  is  zero;  (b)  the  tangential  velocity 
does  not  have  normal  gradient;  and  (c)  the  divergence  of  a  boundary  cell  is 
zero  In  this  research,  an  improved  partial  cell  treatment  is  developed  to 
fulfill  the  above  three  conditions.  The  partial,  cell  is  defined  as  the  cell 
that  intersects  curved  boundary  or  internal  obstacl.e.  The  curved  boundary 
surface  is  approximated  by  a  plane.  The  basic  orientatioxi  of  the  plane  is 
described  by  the  parameter  NFB  through  its  numerical  values  1,  2,  3,  <+,  5,  or 
6,  The  definition  of  numerical  values  for  the  NFB  is  specified  in  Table  3. 

The  value  of  NFB  is  only  dependent  on  mesh  division  and  does  not  change 
tlirough  the  computation.  The  two  slopes  of  the  boundary  plane,  ANB  and  BNB 
can  be  assigned  according  to  the  definition  in  Table  4. 
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TABLE  3 


DEFINITION  OF  VALUES  FOR  NFB 


NFB--1 

obstacle 

:  cell, 

,  all 

.  faces  of  the  cell  are  closed. 

NFB=  1 

partial 

cell , 

the 

right  face  of  the  cell  is  closed. 

NFB-  2 

partial 

cell , 

the 

left  face  of  the  cell  is  closed. 

NFB-  3 

partial 

cell , 

the 

bottom  face,  of  the  cell  is  closed 

NFB-  4 

partial 

cell , 

the 

top  face  of  the  cell  i.s  closed. 

NFB-  5 

partial 

cell , 

the 

front  face  of  the  cell  is  closed. 

NFB-  6 

partial 

cell , 

the 

back  face  of  the  cell  is  closed. 

TABLE  4 

THE 

DEFINITION 

OF  VALUES  FOR 

ANB  AND 

BNB 

NF 

-  1.2 

ANB 

BNB  - 

ax 

dY 

az 

NF 

-  3,4 

ANB-^ 

BNB  - 

ay 

aZ 

ax 

NF 

-  .3,6 

ANB  -  — 

BNB  - 

az 

dX 

ay 

In  the  improved  partial  cell  treatment,  the  normal  and  tangential 
velocity  conditions  on  the  boundary  plane  are  first  considered.  Fig.  7  gives 
an  illustration  for  the  case  when  NFh-2 .  In  this  case,  the  fluid  is  mostly  on 
the  riglit  side  of  the  boundary  plane  and  the  left  side  of  the  boundary  is 
closed.  The  unknown  velocity  component  is  set  so  that  the  zero 
velocity  normal  to  a  boundary  surface  is  fulfilled.  The  normal  velocity 
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can  be  written  as 


Uji  “  ucos(n,x)  +  vcos(n,y)  +  wcos(n,z)  (59) 

Three  direction  cosine,  cos  (n,x),  cos  (n,y)  and  cos(n,z)  are  given  by 


cos(n,x)  - 


ANB^  +  BNB^  +  1 


cos(n,y) 


-MB 


(60) 


MB^  +  BNB^  +  1 


cos(n, z) 


-BNB 


MB^  +  BNB^  +  1 


u,  V,  ai\d  V7  are  the  average  velocities  in  x,  y,  and  z  direction,  respectively, 
and  are  expressed  as; 


(61) 


-  <”i,j,k-i  +  "i,j,k> 

where  A  -  7^/(!ix^;7j^  is  the  distance  between  the  boundary  face  and  the  center 

of  right  face,  as  .shown  in  Fig,  7(a).  By  substituting  Eq .  (60)  and  Eq ,  (61) 

into  Eq .  (59)  and  letting  left  side  of  Eq .  (59)  be  zero,  we  can  obtain 

u.  T  ,  ,  as  follows 
1-1, J  ,k 


u . 


i-l,j  ,k 


1-A 

*i  i  k 
A  2Acos(n,x) 


-u . 


2Acos(n,x) 


■(w. 


.  ,  ,  T  +  w.  .  ,  ;-.:os(n,z) 

i-J-k-l  i,J.k 


(62) 
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Now,  the  continuity  equation  (20)  is  satisfied  by  adjusting  the  partial 
cell  pressure.  The  difference  between  the  partial  cell  and  interior  cell  is 
the  value  of  partial  cell  geometrical  quantities  AR^  j  AT^  j  I,;, 
and  AFi  j  ij.  j  is  the  fractional  volume  of  cell  (i,j,k)  open  to  fluid; 
^i,j,k'  and  j  are  the  fractions  of  the  area  of  the  right  face, 
the  top  face  and  the  front  face  or  cell  (i,j,k)  that  is  open  to  flow.  These 
geometrical  quantities  need  to  calculate  according  to  the  definition. 
Fig.  7(b)  is  the  example  for  the  case  when  NFB“2 .  It  is  noted  that  all  these 
geometrical  quantities  are  simply  set  to  1  for  interior  cell. 

For  the  free  surface  boundary  conditions,  the  stress  tangential  to  the 
surface  must  vanish;  therefore,  the  velocities  of  both  fluids  must  be  equal 
and  the  stress  trormal  to  the  free  surface  must  be  exactly  balanced  by 
externally  applied  normal  stresses. 

In  the  finite  difference  scheme,  the  normal  free  surface  boundary 
condition  is  satisfied  by  setting  the  surface  cell  pressures  Pi,j,k  ®qual  to 
the  value  calculated  by  a  linear  interpolation  between  the  surface  pressure  pg 
and  pressure  pj^  of  the  adjacent  cell.  In  the  absence  of  the  surface  tension, 
the  pressure  at  the  free  surface  is  set  equal  to  the  atmospheric  pressure. 

The  tangential  free  surface  boundary  condition  is  satisfied  by  setting 
velocities  on  every  cell  boundary  between  a  surface  cell  and  an  empty  cell. 

If  the  surface  cell  has  only  one  neighboring  empty  cell,  the  boundary 
velocity  is  set  to  insure  the  vanishing  of  velocity  divergence  defined  in 
Eq.  (20).  VJlien  there  are  two  or  more  empty  neighbor  cells.  the  itidividual 

...  ,  du  dv  dw  , 
contributions  to  tlie  divergence  — ,  — ,  —  are  separately  set  to  zero. 

dx  dy  dz 

Numerical  calculations  often  have  computed  quantities  that  develop  large, 
high-frequency  oscillations  in  space,  time,  or  both.  This  behavior  is  usually 
referred  to  as  a  numerical  instability,  especially  if  die  physical  problem 
being  studied  is  known,  not  to  have  unstable  solutions.  When  the  physical 
problem  does  have  unstable  solutions  and  if  the  calc;ulated  results  exhibit 


significant  variations  over  distances  comparable  to  a  cell  width  or  over  times 
comparable  to  the  time  ii  cremeiit ,  the  accuracy  of  the  results  cannot  be  relied 
on.  To  prevent  this  type  of  numerical  instability  or  inaccuracy,  certain 
restrictions  mvist  be  observed  in  defining  the  mesh  increments  ,  iSyj  and 
the  time  increment  (5t,  and  the  upstream  differencing  parameter  a. 

For  accuracy,  the  mesh  increments  must  be  chosen  small  enough  to  resolve 
the  expected  spatial  variations  in  all  dependent  variables.  The  choice  of  the 
time  increment  necessary  for  stability  is  governed  by  two  restrictions. 
First,  the  convective  limit  (the  Courant  condition).  The  material  cannot  move 
through  more  than  one  cell  in  one  time  step  because  the  difference  equations 
assume  fluxes  only  between  adjacent  cells.  The  numerical  expression  may  be 
written  as 


<5t  <  MIN 


6x . 

1 


<5z 


k 


I V .  ,  I  I  w .  .  , 

'  1 .  J  ,k'  '  1,  J  ,1 


(63) 


where  the  minimum  is  taken  over  all  cells  of  the  computing  mesh,  but  it  is 
usual  to  require  that  6t  be  no  more  than  a  small  fraction,  e.g,,  1/4  of  the 
minimum  cell  transit  time. 

Second,  the  diffusive  limit  states  that  when  a  non-zero  value  of 
kinematic  viscosity  is  used,  momentum  must  not  diffuse  more  than  approximately 
one  cell  in  one  time  step.  A  linear  stability  analysis  shows  that  this 
limitation  implies 


vdt  <  —  ( - +  -  -I-  -  )  (64) 

2  ^2  ,2  _  2 

(5x.  oy .  oz. 

L  j  k 

with  6t  chosen  to  satisfy  the  above  two  inequalities,  the  parameter  q  describ¬ 
ing  the  proportion  of  donor  cell  differencing  should  have 


1  >  a  >  MAX  i 


' 

A 

u.  .  ,  dit 
1  .  .  k 

V .  .  ,  (5t 

t ,  J  .  k 

w.  .  , 6t 
i.J  .k 

1 

(5x . 

1 

J 

t 

<5z, 

k 

((^5) 
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As  a  rule  of  thumb,  an  a  approximately  1.2  to  1.5  times  larger  than  the  right' 
hand  member  of  the  last  inequality  is  a  good  choice.  If  a  is  too  large,  an 
unnecessary  amount  of  numerical  smoothing  (diffusion-like  truncation  errors) 
may  be  introduced. 

This  improvement  of  the  volume  of  fluid  (VOF)  technique  allows  the 
numerical  simulation  of  three-dimensional  liquid  sloshing  in  a  container  of 
arbitrarj'  geometry.  Major  improvements  were  the  taking  into  account  of  free 
surface  orientation,  transporting  hexahedral  shape  fluid  volume  from  cell  to 
cell  and  considering  the  normal  and  tangential  velocity  boundary  conditions  on 
curved  solid  boundaries  in  partial  cell  treatments.  The  following  chapters 
describe  simulation  results  in  various  application  conditions. 
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CHAPTER  3;  SUBHARMONIC  RESPONSES  UNDER  VERTICAL  EXCITATION 
AND  ITS  CONTROL  BY  BAFFLESINTRODUCTION 

It  is  known  that  vertical  vibration  causes  a  quiescent  liquid  surface  to 
become  unstable  and  the  frequency  of  surface  oscillations  usually  occurs  at 
exactly  one-half  that  of  the  container  motion^^^) .  As  part  of  the  algorithm 
testing,  numerical  experiments  are  carried  out  to  examine  certain  phenomena 
associated  with  flow  instability. 

A  rectangular  tank  of  25.4-  cm  long,  25.4  cm  wide  and  35.56  cm  high  is 
used  for  the  numerical  simulation.  Two  different  water  depths(^)  are  chosen, 
h  -  19.05  cm  for  deep  water  care  and  h  -  6.35  cm  for  shallow  water  sloshing. 
Zero  initial  velocity  field  and  hydrostatic  pressure  distribution  are  fir.st 
assumed.  Later,  a  small  velocity  perturbation  with  horizontal  component 
-  0.127  cm/s  and  vertical  component  Vq  -  0.127  cm/s  are  introduced  to 
accelerate  the  growth  of  instability.  In  all  cases,  vertical  harmonic 
excitations  yoSin  ut  are  applied  v;ith  maximum  amplitude  of  y^  1,27  cm  and 
the  exciting  frequency  slightly  greater  than  twice  of  the  first  linear 
natural  sloshing  frequency 
Deep  Water  Sloshing 

Sloshing  in  the  tank  without  baffles.  The  frequency  of  forced  oscillation 
is  chosen  as  22.258  rad/s  (2.04coj^).  The  computed  free  surfaces  are  shown  in 
Fig.  8(a)  through  Fig.  8(e)  at  five  different  times;  2  . 9  5  s  (,  10 .  5T )  , 
3 . 24s(ll . 5T) ,  3,81s(13.5T)  ,  4,09s(14.5T)  and  4 . 3 7s ( 15 . 5T) .  T  is  the  period  of 
forced  oscillation.  Fig.  8(f)  shows  the  free  surface  at  4,37s(15,5T)  without 
higher  order  initial  velocity  perturbation.  The  computed  velocity- vector 
field  on  vertical  (x,y)  and  (z,y)  planes  are  shown  in  Fig.  9.  In  the  figures. 
Section  I  and  K  are  parallel  to  the  left  wall  of  the  tank  and  parallel  to  the 
front  wall  of  the  tank  respectively.  The  subharmonic  responses  of  free 
surface  are  indeed  observed  in  the  simulation  results.  That  is,  the  surface 
wave  only  oscillates  half  period  when  the  forced- oscillation  passes  one 
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period.  Since  the  iiatural  frequency  of  free  sloshing  in  x  direction  is  equal 
to  that  in  z  direction  due  to  the  tank  configuration,  the  subharraonic 
responses  occur  in  the  both  directions  at  the  same  time.  As  a  result,  the 
free  surface  waves  are  obviously  three-dimensional,  with  maximum  wave  height 
on  the  tank  corner.  It  is  noted  that  vertical  oscillation  causes  a  quiescent 
free  surface  to  become  unstable.  Once  this  occurs,  the  amplitudes  of  surface 
wave  increase  rapidly.  In  this  numerical  example,  the  surface  wave  finally 
impacts  the  tank  top  and  then  breaks  with  bubbles  in  the  liquid. 

The  higher-order  initial  velocity  perturbation  is  used  to  stimulate  the 
dynamic  instability  of  vertical  sloshing.  The  numerical  results  proved  this 
to  be  very  effective.  Ti\e  idea  of  using  higher-order  perturbation  is  also 
reasonable  because  there  always  exists  perturbations  in  physical  problems. 

SloshinK  in  the  tank  with  horizontal  splitter  ring.  The  horizontal 
splitter  ring  is  placed  at  the  distance  of  14  cm  above  the  tank  bottom.  In 
numerical  simulation,  the  splitter  ring  is  created  by  blocking  out  the  thirty- 
six  appropriate  computational  cells.  The  computed  free  surfaces  are  shown  in 
Fig.  10  at  two  different  times;  2.95s(10.5T)  and  4 . 37s(15 . 5T) .  Fig.  11  shows 
the  computed  velocity-vector  field  on  vertical  (x,y)  and  (z,y)  planes.  In  the 
numerical  results,  only  small  wave  motion  can  be  observed.  This  indicates 
that  surface  responses  have  been  suppressed  by  splitter  ring. 

Sloshing  in  the  tank  with  vertical  splitter _ ring.  The  vertical  splitter 

ring  is  so  placed  that  it  is  parallel  to  the  front  wall  of  the  tank  and  11.43 
ciii  fruiii  the  front  wall.  In  numerical  siiuulation,  the  splitter  ring  is  created 
by  blocking  out  the  forty-six  appropriate  computational  cells.  The  computed 
free  surfaces  are  shown  in  Fig,  12  at  three  different  times;  2 , 95s ( 10 . 5T)  , 
3.24s(11.5T)  and  4 . 37s ( 15 . 5T) .  Fig.  13  shows  the  computed  velocity- vector 
field  on  vertical  (x,y)  and  (z,y)  planes.  The  computed  free  surfaces  still 
oscillate,  but  tlie  responses  of  free  surface  decrease  largely  in  the  direction 
perpendicular  to  tlie  splitter  ring.  So,  the  free  surface  wave  oscillates 
mainly  in  the  direction  parallel  to  the  splitter  ring. 
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Shallow  Water  Sloshing 


Sloshing  in  the  tank  without  baffles.  The  frequency  of  forced  oscillation 
is  chosen  as  18.188rad/s  (2.04w^).  The  computed  free  surfaces  are  shown  in 
Fig.  14(a)  through  Fig.  14(e)  at  five  different  times;  5 . 69s ( 16 . 5T) , 
6.04s(17.5T) ,  6. 74(19. 5T),  7.08s(20.5T)  and  7 .42s(21 . 5T) .  Fig.  14(f)  shows 
the  free  surface  response  at  7.42s(21.5T)  without  higher-order  initial 
perturbation.  The  computed  velocity-vector  field  on  vertical  (x,y)  and  (z,y) 
planes  are  shown  in  Fig.  15  through  17.  From  the  numerical  results,  the 
subharmonic  responses  also  occur  in  shallov;  water  depth.  In  addition,  with 
large  surface  motion,  some  parts  of  tank  bottom  become  uncovered. 

Sloshing  in  the  tank  with  horizontal  splitter  ring.  The  location  of 
horizontal  ring  is  the  same  as  in  the  deep  v;ater  case.  The  computed  free 
surfaces  are  shown  in  Fig.  18  at  four  different  times;  5 . 6 9s ( 16 . 5T) , 

6 . 04s ( 17 . 5T)  ,  7.08s(20,5T)  and  7 . 42s (21 . 5T) .  The  computed  velocity-vector 
field  on  vertical  (x,y)  and  (z,y)  planes  are  shown  in  Fig.  19  through  Fig.  21. 
The  surface  responses  are  the  same  as  with  the  no  baffle  case  as  long  as  the 
surface  wave  does  not  impact  the  splitter  ring.  The  surface  wave  does  not  go 
beyond  the  splitter  ring,  even  chough  the  large  surface  wave  impacts  it. 

Sloshing  in  the  tank  with  vertical  splitter  ring.  The  location  of 
vertical  splitter  ring  is  the  same  as  in  the  deep  water  case.  The  computed 
free  surfaces  are  shown  in  Fig.  22  at  three  different  times;  5 . 69s ( 16 . 5T) , 
6 , 04s ( 17 . 5T ) ,  and  7 , 42 s ( 2 1 . 5T )  .  The  computed  velocity- vector  field  on 
vertical  (x,y)  and  (z,y)  planes  arc  shown  in  Fig.  21.  In  this  case,  the  part 
of  splitter  ring  on  the  tank  bottom  has  much  influence  on  free  surface  motion. 
Therefore,  the  free  surface  response  on  the  direction  perpendicular  to 
splitter  ring  is  much  less  in  comparison  with  deep  water  case. 
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CHAPTER  4;  THE  CONTROL  OF  LATERAL  SLOSHING  IN  RECTANGULAR 
TANKS  -  THE  EFFECTIVENESS  OF  FIXED  BAFFLE  AND 
MOVING  BAFFLE  WITH  FEEDBACK  CONTROL 

Within  the  moving  container,  various  types  of  sloshing  waves  can  be 
created  depending  on  the  liquid  depth  and  motion  frequeiicy.  Sloshing  is 
essentially  a  nonlinear  phenomena.  For  large  amplitude,  excitations,  the 
liquid  response  can  be  rather  violent.  Baffles  have  been  used  to  suppress  and 
control  sloshing.  However,  it  may  be  desirable  to  explore  the  possibilities 
of  using  active  moving  baffles  to  control  sloshing  when  the  tank  acceleration 
is  large.  The  subject  study  was  to  investigate,  numerically,  the  relative 
effectiveness  of  passive  (fixed)  baffles  and  active  (moving)  baffles  for 
control  of  liquid  sloshing.  Both  the  deep  water  case  and  the  shallow  water 
case  were  studied.  The  study  showed  that  active  control  of  liquid  sloshing, 
using  moving  baffles  with  feedback  mechanisms,  can  be  very  effective  in 
suppressing  large  amplitude  sloshing. 

Fig.  24  illustrates  the  concept  of  active  baffle  for  the  sloshing 
control.  The  baffle  is  set  up  to  move  in  the  opposite  sense  of  certain 
neighboring  flow  to  suppress  the  sloshing.  In  the  deep  water  case,  a  vertical 
component  of  the  baffle  movement  is  provided  while  the  horizontal  baffle 
velocity  is  given  for  the  shallov/  water  case.  Issues  in  numerical 
implementation  included  the  appropriate  set-up  of  boundary  condition  for 
moving  baffles,  chocsing  the  appropriate  coiitrol  strategy,  such  as  using 
suitable  information  to  provide  feedback  and  choosing  appropriate  feedback 
coefficient.  Furthermore,  noisy  data  needs  to  be  smoothed  out  to  provide 
feedback . 

For  the  horizontal  baffles  in  Fig.  24(a) ,  the  feedback  was  provided  by 
V*,  the  vertical  velocity  of  the  cell  near  the  baffle,  For  the  left  baffle  at 
(lOBS,  JOBS,  KOBS) ,  the  baffle's  vertical  velocity  is 

vb  =  -  .  V-  [  ^  1 

2 


(66) 
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where  v  is  the  coefficient  of  feedback.  For  the  vertical  baffle  in 
Fig.  24(b),  the  feedback  was  controlled  through  U*,  the  horizontal  velocity  of 
the  cell  near  the  baffle,  the  baffle's  horizontal  velocity  is 


Ub 


V 


Y(J0BS)-Y(2) 

1  A  y 

2 


+  1 


(67) 


Since  the  computed  time  history  of  the  reference  velocity  exhibits  noisy 
character  of  complicated  three-dimensional  nonlinear  wave,  cross-validation 
method  is  used  to  e.stimate  the  saioothing  parameter.  The  treated  reference 
velocity  information  was  extrapolated  to  provide  the  feedback  for  active 
baffle.  The  coefficient  of  feedback  u  is  defined  in  Fig.  24.  The  case  of 
1^=0,  represents  fixed  baffle. 

A  rectangular  tai^k  of  38.1  cm  long,  38.1  cm  wide,  and  38.1  cm  high  is 
u.sed  for  the  nujuerical  simulation.  Water  depth  h-28.6  cm  is  chosen  for  deep 
water  sloshing  simulation.  Two  baffles,  one  on  each  side,  are  located  near 
equilibrium  free  surface.  Horizontal  harmonic  excitations  xo'“3.048  cm  and  the 
excitation  frequency  w  equal  to  1.01  times  the  first  linear  natural  sloshing 
frequency.  Fig.  25(a-b)  shows  the  free  surface  plot  at  t“8.25T  for  the  case 
when  the  baffles  are  fixed  and  the  case  when  the  baffles  are  active  with 
feedback  coefficient  j^-0.3  after  a  first  quarter  period.  The  wave  height  time 
history  on  the  left  wall  of  each  tank  is  shown  in  Fig.  26.  The  scale  is  given 
in  inches,  and  the  tank  top  at  Sa’=15  is  indicated  as  a  dashed  line.  The 
effect  of  active  baffle  is  obvious  wave  height  reduction  of  50%  was  indicated. 
The  shallow  water  computation  is  carried  out  at  a  25%  fill-depth  for  the  same 
tank  with  same  amplitude  for  a  forcing  frequency  at  1.05  times  the  fundamental 
linear  natural  frequency.  A  single  baffle  is  located  at  the  center. 
Fig.  27(a-b)  shows  the  free  surface  plot  for  fixed  and  active  baffle  (i/“0.5), 
respectively.  Drastic  reduction  of  wave  height  on  the  left  wall  (more  than 
60%),  after  the  activation  of  moving  baffle  at  t=T/4 ,  can  be  seen  from 
Fig.  28 (a-b) . 
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CHAPTER  5;  THE  IMPACT  IN  1  g  ENVIRONMENT  BY  BAFFLESINTRODUCTION 


For  the  study  of  impact  in  1  g  environment,  breaking  wave  against 
vertical  wall  seated  on  a  structural  foundation,  was  numerically  simulated. 
Laitone's  second  order  solitary  wave  theory  is  used  as  the  initial  conai tions . 
When  no  foundation  is  present,  the  rvin  up  of  a  solitary  wave  on  a  vertical 
wall  for  a  range  of  wave  height/water  depth  ratio  agrees  well  with  the 

t 

experimental  data.  The  phenomena  of  wave  breaking,  such  as  wave  steepening, 
overturning  and  formation  of  bores  have  been  successfully  simulated.  Very 
high  intensity  shock  pressure  and  wave  impact  force  on  the  vertical  wall  are 
also  obtained. 

A  mesh  of  180  cells  in  the  x-direction  and  50  cells  in  the  y  direction 
was  used  to  represent  the  computation  region.  The  space  increments  (Sx^Acm  and 
(Sy— 2cm  in  x  and  y  direction  were  used  for  all  cases,  Definition  sketch  is 
shown  in  Fig.  29.  Water  depth  d  is  35cm;  initial  wave  cre.st  i.s  located  at 
Xo“240cra:  water  density  p-lg/cm^ ;  gravity  acceleration  g-980cm/s2;  kinematic 
viscosity  i/“0 . 01002cm^/s  ;  iteration  convergence  cri.teri.on  f-0.001. 

First,  solitary  wave  propagating  toward  a  vertical  wall  without  formation 
is  computed.  The  wave  run-up  ratio  R/d  is  compared  with  the  experimental  data 
of  Street  &  Camfield  (1988)  [40]  in  Fig.  30.  Clearly,  the  numerical  results 
are  in  excellent  agreement  with  experiments. 

The  computations  have  been  carried  out  for  several  configurations  of 
foundations.  The  berm  of  the  foundation  is  160cm  and  slope  is  1:2.  The 
relative  water  depth  at  the  berm  of  foundations  d]^/d  are  0.^1,  0.60,  0.49, 
0,37  and  0.2.  Wave  height  H  is  18 . 9cm(H/d=0 . 54)  .  Flow  field  for  shallow 
foundation  (dp/d=0.71)  is  shown  in  Fig,  31.  In  this  case,  the  solitary  wave 
is  fully  reflected  by  the  wall  without  breaking.  Fig.  32  is  th-  flow  field 
for  one  of  middle  size  foundation  (d]/d=0.37).  As  a  wave  propagates  on  the 
berm,  water  particle  velocity  at  the  wave  crest  increases  and  wave  front 
becomes  steeper.  It  eventually  becomes  unstable  and  breaks  when  water 
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particle  velocity  at  the  wave  crest  exceeds  the  wave  speed  (theoretical  value 
is  230cni/s),  Flow  field  for  fairly  high  foundation  (dj_/d“0.2)  is  showit  in 
Fig.  33.  It  shows  clearly  the  plunging  breaker  in  front  of  the  wall.  After 
breaking,  the  bores  are  formed  with  high  speed.  The  shock  pressure  occurs 
when  the  botes  hit  the  wall.  Fig.  34  and  35  are  the  evolution  of  waves  for 
di/d“0.37  and  0.2.  The  time  history  of  pressure  at  still  water  level  on  the 
wall  is  shown  in  Fig.  36.  From  Table  5,  it  is  easier  to  find  the  maximum 
shock  pressure  and  impulse  force.  Shock  pressure  on  the  wall,  in  the  .  ase  of 
middle  site  foundatiox-i,  is  smaller  than  that  in  the  case  of  high  foundat,.on, 

TABLE  5 

MAXIMUM  WAVE  PRESSURE  AND  FORCE  VALUES  (H/d«0.54) 


dl/d 

1.0 

0.71 

0.60 

0.49 

0.37 

0.2 

^'max/ 

1.48 

1 . 54 

1.73 

1 . 99 

2 . 57 

3 , 31 

^max/ 

2.14 

2.75 

3.15 

4.01 

5.36 

15.51 

Fig.  3  7  to  Fig.  39  show  the  comput.ational  results  for  wave  height 
H”27 . 3cm(H/d"0 , 7 8 )  whicti  is  the  limiting  wave  height  from  solitary  wave 
theory.  Again,  it  shows  the  plunging  breaker  in  front  of  wall.  The  large  jet 
is  ejected  forward  from  the  tip  of  the  wave.  Very  high  shock  pressure  occurs 
v;hen  bore  front  hits  tiie  wall.  The  .shock  pressure  5.32  and 


impulse  force  F,f_^.jj^/pgUd  is  8,/0. 


CHAPTER  6:  SWIRLING  IN  CYLINDRICAL  TANK  ftUD  ITS  CONTROL 


It  is  well  knovm  that  lateral  excitation,  near  the  lowest  liquid  resonant 
frequency,  causes  an  interesting  type  of  liquid  instability  to  take  place  in 
the  form  of  swirl  motion  superimposed  on  the  normal  sloshing  motior  As 
described  in  Reference  1,  "The  motion  is  even  more  complicated  as  a  ty 
'beating'  also  seems  to  exist;  the  first  antisynunetric  liquid- sloshing  mot  . 
first  begins  to  transform  itself  into  a  rotational  motion  increasing  in 
angular  velocity  in,  say,  the  counterclockwise  direction,  which  reaches  a 
maximum  and  then  decreases  essentially  to  zero  and  then  reverses  and  increases 
in  angular  velocity  in  the  clockwise  direction,  and  so  on,  alternately .  The 
frequeitcy  of  rotation  is  less  than  that  of  the  surface  wave  motion  and, 
therefore ,  the  liquid  appears  to  undergo  a  vertical  up-and-down  motion  as  it 
rotates  about  the  tank  axis;  the  rotational  frequency  about  this  up-and-down 
axis  is  about  the  same  as  that  or  the  wave  motion."  As  part  of  the  algorithm 
testing,  numerical  experiments  were  carried  out  to  examine  this  phenomena  of 
flow  instability. 

An  upright  cylindrical  tank  of  diameter  d-19,558  cm  and  water  depth  h-d, 
is  used  for  the  numerical  simulation.  The  dimensions  of  tlie  cylindrical  tank 
for  numerical  simulation  are  the  same  as  the  experimental  tank  in  Ref.  91. 
With  this  diameter  and  water  deptli,  the  first  ant  isyinme  tr  ical  natural 
frequency  Wi("13.574  rad/s  (u?!d/g"3 . 68)  .  The  amplitude  of  the  forced- lateral 
oscillation  is  Xo=0.3363  cin(xc/d=0 . 0172)  .  The  frequency  of  forced  oscillation 
i..s  W'"13.979  rad/s  (w2d/g=‘3 , 90)  .  The  parameters  of  f  or  ced  -  osc  i  1  lat  ion  are 
selected  in  the  sv/irl  region  predicted  by  lahoratory  experiment  (Fig.  40),  so 
that  the  swirl  motion  of  liquid  c.an  be  studied  by  numerical  simulation. 

The  liquid  is  initially  at  rest.  Under  a  sinusoidal  excitation,  the  time 
history  of  wave  lieight,  pressure  and  force  components  are  shown  in  Fig.  41. 
The  location  for  v/ave  height  record  Sg  is  4.191  cm  from  the  left  tank  v/all. . 
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The  location  for  pressure  record  is  4.191  cm  from  the  left  tank  wall  and 
18.161  cm  from  the  tank  bottom.  force  component  in  the  direction  of 
excitation  and  is  force  component  perpendicular  to  the  paper  plane 
(transverse  force).  "Beai'ing"  is  apparent.  Fig.  42  shows  computed  free 
surface  within  one  period  of  forced  oscillation  T;  the  time  interval  is 
0.125  T  The  computed  velocity-vector  field  on  horizontal  (x,z)  plane,  is 
shown  in  Fig.  43,  in  which  section  J"S  and  J“9  are  18.161  cm  and  20.955  cm 
from  the  tank  bottom,  respectively.  The  computed  velocity-vector  field  on  Che 
vertical  (x,y)  and  (z,y)  plane  is  shown  in  Fig.  44.  From  the  numerical 
results,  swirl  motion  of  free  surface  can  be  observed,  in  agreement  with 
experiment  results  of  Ref.  41.  The  period  of  swirling  motion  i.s  the  same  as 
the  period  of  forced  oscillation.  The  horizojital  velocity -vector  field  shows 
the  swirling  direction.  One  important  characteristic  of  liquid  sX'/irl  motion 
i.s  that  large  transverse  force  component  exists.  This  may  be  used  to 
determine  whether  or  not  the  swirl  motion  of  free  surface  occur.s . 

Systematic  n\JLmerical  experiment.^  were  carried  out  to  compare  numerical 
predicted  wave  amplitude  with  laboratory  test  results  of  Ref.  41.  Four  levels 
of  maximum  aniplicude.s  of  o.scillanion,  Xo/d=0.0056,  0.0115,  0.0172  and  0.0227, 
were  used.  In  the  numerical  model,  tbie  number  of  cells  in  x,  y  and  z 
direction  are  7,  12  and  7  respectively.  The  position  for  wave  height  record 
is  1.397  ciTi  from  the  left  tank  wall.  Liquid  amplitude  re.spon.ses  in  the  first 
aritisymiiietr ic  slosh  mode  ate  .shown  a.s  a  function  of  excitation  frequency,  for 
several  constant  values  of  excitation  amplitude  Xo/d  in  Fig.  40,  For  eaci, 
excitation,  the  numerical  data  is  obtained  by  averaging  successively  10  v/ave 
heights  after  steady -state  .solutioii  is  approached,  The  numerical  result-s 
agree  well  with  experimental  ilata  when  ex.citation  frequency  i.s  less  than  0^ , 
When  excitation  frequency  is  larger  chan  Wj/,  the  numerical  re.sults  predict  a 
slightly  larger  result,  but  it  .still  shows  good  agreement  with  experimental 
results.  In  the  shaded  region,  the  obviously  sv/irl  motion  occurs  as  jjredicted 
by  no.iierical  simulation. 
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EFFECTS  OF  BAFFLES  IN  THE  LATERAL  SLOSHING  OF  A  CYLINDRICAL  TAiif; 


The  introduction  of  baffle  changes  the  characters  of  liquid  sloshing  in 
many  complicated  ways.  To  evaluate  baffle  effect,  four  typical  baffle 
arrangements  were  each  added  to  the  baseline  case  described  earlier  (Figs.  41- 
44)  and  tested  numerically.  Shown  in  Fig.  45,  is  the  arrangement  of  these 
baffle  types,  namely, 

(a)  Short  vertical  splitter  plates  in  transverse  direction. 

(b)  Long  vertical  splitter  plates  in  transverse  direction. 

(c)  Short  vertical  splitter  plates  in  the  direction  of  excitation. 

(d)  Horizontal  ring  at  undisturbed  free  surface. 

Short  Vertical  Splitter  Plates  in  Transverse  Direction 

The  arrangement  is  shown  in  Fig.  45(a).  The  time  history  of  wave  height, 
pressure  and  force  components  is  shown  in  Fig.  46.  Fig.  47  shows  the  computed 
free  surface  within  one  period  of  forced  osci.llation .  The  computed  velocity- 
vector  field  on  horizontal  (x,z)  plane,  is  shown  in  Fig.  48.  The  computed 

velocity-vector  field  on  vertical  (x,y)  and  (z,y)  plane  is  shown  in  Fig.  49. 
Swirl  motion  of  free  surfaces  are  observed  in  this  configuration  and  there 
still  exists  a  large  transverse  force  shown  in  Fig.  48(d).  However,  the 

amplitude  of  free  surface  wave  is  reduced  due  to  the  effect  of  splitter 
plates  . 

Long  Vertical  Splitter  Plates  in  Transverse  Direction 

The  arrangement  is  showti  in  Fig.  45(b).  The  time  history  of  wave  height, 
pressure  and  force  component.s  is  shown  in  Fig.  50.  Fig.  51  shows  the  computed 
free  surface  within  one  period  of  forced  o.scillation .  The  computed  velocity- 
vectui.  fielci  on  horizontal  (x,z)  plane  i.s  shown  in  Fig,  52.  The  computed 

veloc  i  ty -vec  tor  field  on  the  vertical  (x,y)  and  (z,y)  plane  i.s  shown  in 
Fig.  53,  Swirl  motion  of  free  surface  is  suppressed  in  this  configuration. 
Only  small  surface  V7ave  can  be  observed.  Small  transverse  force  .shown  in 
Fig.  50(d)  indicates  no-.swirl  motion. 
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Short  Vertical  Splitter  Plates  in  the  Direction  of  ExciV-ation 


The  baffle  type  is  shown  in  Fig.  45(c).  The  time  history  of  wave  height, 
pressure  and  force  components  is  shown  in  Fig.  54.  Fig.  55  shows  the  compmted 
free  surface  within  one  period  of  forced  oscillation.  The  coioputed  velocity- 
vector  field  on  horizontal  (x,z)  plane  are  shottjn  in  Fig.  56.  The  computed 
velocity-vector  field  on  the  vertical  (x.yO  and  (z,y)  plane  is  shown  in 
Fig.  57.  Although  swirl  motion  of  free  surface  is  easily  suppressed  in  this 
configuration,  very  large  wave  motion  in  the  direction  of  excitation  is 
formed.  Only  small  transverse  force  is  shown  in  Fig.  54(d). 

Horizontal  Rina  at  undisturbed  Free  Surface 

The  ring  is  located  as  shown  in  Fig.  45(d).  The  time  history  of  height, 
pressure  and  force  components  is  shown  in  Fig.  58.  Fig.  59  shoves  the  computed 
free  surface  within  one  period  of  forced  oscillation.  The  computed  velocity- 
vector  field  on  horizontal  (x,z)  plane  is  shown  in  Fig.  60.  The  computed 
velocity- vector  field  on  the  vertical  (x,y)  and  (z.y)  plane  is  shown  in 
Fig.  61,  Swirl  motion  of  free  surface  is  not  obviou.s  in  this  configuration. 
But,  one  can  conclude  that  there  still  exists  swirl  motion  from  large 
transverse  force,  as  shown  in  Fig.  58. 

Swirling  and  Drainage 

Various  combinations  of  swirling  and  drainage  v/ere  computed.  The 
swirling  intensification  and  the  air - entra inmenc  induced  during  the  liquid 
drainage  through  the  tank  bottom  v/ere  noted.  Surface  and  surface  contour 
plots  at  various  rates  of  drainage  were  shown  in  Fig.  62.  A  preliminary 
laboratory  experiment  was  carried  out  which  revealed  rather  complicated 
instability  proce.sses  involved  in  the  drainage.  To  include  such  processes 
would  require  further-refined  numerical  and  laboratory  modeling.  Therefore, 
the  topic  was  only  briefly  mentioned  in  this  report. 
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CHAPTER  7:  CONCLUSION  AND  FUTURE  WORK 


Fuel  sloshing  in  space  environment  is  a  problem  of  technica]. 
significance.  The  low  gravity  environment,  the  rapid  maneuvering  of  the 
spacecraft,  the  flexible  container,  impact  loads,  surface  tension,  flow- 
induced  sloshing  (such  as  drainage)  arid  the  sloshing  control  and  suppression 
offer  many  areas  of  possibility  of  advancement.  Some  progresses  in  sloshing 
simultion  and  control  were  reported  in  the  present  report.  Many  problems 
remain,  relating  to: 

•  Low  gravity  environment 

•  Surface  tension 

•  Sub -grid  phenomena 

•  Small  scale  flow  physics  and  flow  instability 

•  Impact  on  flexible  structure 

•  Improving  surface  representation 

•  Experimental  verification 

More  studies  are  required  to  address  these  issues. 
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Fig.  10.  Surface  plot  for  deep  water  sloshing  under  vertical  excitation 
(horizontal  baffles) 
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Fig.  12.  Surface  plot  for  problem  of  deep  water  sloshing  under  vertical 
excitation  (vertical  baffles) 
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Fig.  I'i.  Velocity  plot  for  problem  of  deep  water  .slo:;Viing  under  vertical 
excitation  (vertical  baffles) 


M 


iP«rturb«tloA 


SECTION  K  •  2 


SECTION  K  ■  3 


section  K  ■  « 


SECTION  K  ■  II 


SECTim  K  • 


Fig.  20,  Velocity  plot  for  shallow  water  sloshing  under  vertical  excitation 
(horizontal  battles) 
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Fig.  21.  Velocity  plot  tor  shallot;  water  sloshing  under  vertical  excitati.on 
(horizontal  baffle.s) 
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Fig.  26.  Effect  of  baffles  on  surface  height  history 
(a)  baffle  fixed;  (b)  baffle  active 


Fig.  27. 


Free  surface  plot  for  shallow  water  cases 
(a)  baffle  fixed;  (b)  baffle  active 


Fig.  28.  Effect  on  baffles  on  surface  height  history 
(a)  baffle  fixed;  (b)  baffle  active 


Fig.  29.  Definition  sketch 
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Fig.  40.  Comparison  of  the  liquid  free  surface 

response,  numarical  results,  - - 

experimental  results  (Ref.  41) 
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Fig.  41.  Wave  height,  pressure  and  force  history  •  lateral  sloshing  in 
cylindrical  tank 


Fig.  51.  Surface  plo!;  showing  the  effect  of  longer  vertical  splitter  plates 
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Fig.  52,  Flow  field  plot  for  horizontal  sections  J«8  and  J-9 
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Fig.  53,  Flow  field  plot  for  vertical  sections  K“5  and  1-5 
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Fig.  55.  Surface  plot  shov;iug  the  effect  of  effect  of  in-line  vertical 
splitter  plates 
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